﻿using System;
using System.Collections.Generic;
using System.Linq;

namespace PostProcCalc
{

    public class MeshFEM : IMeshFEM
    {

        #region Private fields

        private List<Point> _points;
        private List<int[]> _nvtr;
        private double[] _q;
        private String[] _matNames;

        private BasisType _basis;
        private ElementType _elements;

        #endregion

        #region Implementation of IMeshFEM

        public ValueType CurrentValueType
        {
            get;
            set;
        }

        public double GetValueInPoint(Point p)
        {
            for(var i=0; i<_nvtr.Count; ++i)
            {
                if(!PointInsideElement(p, i, _elements))
                    continue;
                return GetValueInPoint(p, i,_elements);
            }

            throw new Exception("Point doesn't belong to area");
        }        

        public IEnumerable<double> GetValuesInPoints(IEnumerable<Point> points)
        {
            return points.Select(GetValueInPoint);
        }

        public double GetDeriveValueInPoint(Point p, ValueType type)
        {

            if (type == ValueType.XDerive || type == ValueType.YDerive)
            {
                for (var i = 0; i < _nvtr.Count; ++i)
                {
                    if (!PointInsideElement(p, i, _elements))
                        continue;
                    return GetDeriveValueInPoint (p, i, type, _elements);
                }

                throw new Exception("Point doesn't belong to area");
            }

            throw new Exception("type is not derivative. ");
        }

        public void Load(params string[] files)
        {
            var fileWorker = new FileWorker(files[0]);

            _points = fileWorker.LoadPoints(files[1]).ToList();
            _nvtr = fileWorker.LoadElements(files[2]).ToList();
            _basis = (BasisType)fileWorker.BasisCode;
            _elements = (ElementType)(_nvtr.First().Count() - 1);
            _q = fileWorker.LoadSolution(files[3]);
            _matNames = fileWorker.LoadMaterialNames(files[4]);

            if (_elements == ElementType.Rectangle)
                RebuildMesh();
            else
                FragmentMesh();
        }

        #endregion

        private double GetValueInPoint(Point p, int numEl, ElementType type)
        {
            double value = 0;
            if (type == ElementType.Rectangle)
            {
                value = GetValueInPointInRectangle(p, numEl);
            }
            else if (type == ElementType.Triangle)
            {
                value = GetValueInPointInTriangle(p, numEl);
            }
            
            return value;
        }

        private double GetDeriveValueInPoint(Point p, int numEl, ValueType valType, ElementType elemType)
        {
            if (elemType == ElementType.Rectangle)
            {
                return GetDeriveValueInPointRectangle (p, numEl, valType);
            }
            else if (elemType == ElementType.Triangle)
            {
                return GetDeriveValueInPointTriangle (p, numEl, valType);
            }
            else
                throw new Exception("value type is not derivative. ");
        }

        private double GetDeriveValueInPointTriangle(Point p, int NumEl, ValueType valType)
        {
            if (valType != ValueType.Value && valType != ValueType.ScalarDerive)
            {
                double detD;

                int[] _elem;

                _elem = _nvtr[NumEl];
                double a01, a02, a03, a11, a12, a13, a21, a22, a23;

                detD = (_points[_elem[1]].X - _points[_elem[0]].X) * (_points[_elem[2]].Y - _points[_elem[0]].Y) - (_points[_elem[2]].X - _points[_elem[0]].X) * (_points[_elem[1]].Y - _points[_elem[0]].Y);

                a11 = (_points[_elem[1]].Y - _points[_elem[2]].Y) / detD;
                a12 = (_points[_elem[2]].Y - _points[_elem[0]].Y) / detD;
                a13 = (_points[_elem[0]].Y - _points[_elem[1]].Y) / detD;

                a21 = (_points[_elem[2]].X - _points[_elem[1]].X) / detD;
                a22 = (_points[_elem[0]].X - _points[_elem[2]].X) / detD;
                a23 = (_points[_elem[1]].X - _points[_elem[0]].X) / detD;

                double derive = 0;

                if (valType == ValueType.XDerive)
                {
                    return _q[_elem[0]] * (a11) +
                           _q[_elem[1]] * (a12) +
                           _q[_elem[2]] * (a13);
                }
                else
                {
                    return _q[_elem[0]] * (a21) +
                           _q[_elem[1]] * (a22) +
                           _q[_elem[2]] * (a23);
                }
            }

            throw new Exception("value type is no derive. ");
        }

        private double GetDeriveValueInPointRectangle(Point p, int NumEl, ValueType valType)
        {
            if (valType != ValueType.Value && valType != ValueType.ScalarDerive)
            {
                int[] _elem;
                _elem = _nvtr[NumEl];

                Point
                    p0 = _points[_elem[0]],
                    p1 = _points[_elem[1]],
                    p2 = _points[_elem[2]];

                double hx = p1.X - p0.X;
                double hy = p2.Y - p0.Y;                

                if (valType == ValueType.XDerive)
                {
                    return _q[_elem[0]] * (1 / hx) * ((p.Y - p2.Y) / hy) +
                        _q[_elem[1]] * (1 / hx) * ((p2.Y - p.Y) / hy) +
                        _q[_elem[2]] * (1 / hx) * ((p0.Y - p.Y) / hy) +
                        _q[_elem[3]] * (1 / hx) * ((p.Y - p0.Y) / hy);
                }
                else
                {
                    return _q[_elem[0]] * ((p.X - p1.X) / hx) * (1 / hy) +
                        _q[_elem[1]] * ((p0.X - p.X) / hx) * (1 / hy) +
                        _q[_elem[2]] * ((p1.X - p.X) / hx) * (1 / hy) +
                        _q[_elem[3]] * ((p.X - p0.X) / hx) * (1 / hy);
                }
            }

            throw new Exception("value type is no derive. ");
        }

        private bool PointInsideElement(Point p, int numEl, ElementType type)
        {
            bool val = false;
            if (type == ElementType.Rectangle)
            {
                val = PointInsideTriangle(p, numEl);
            }
            else
            {
                val = PointInsideRectangle(p, numEl);
            }

            return val;
        }

        private bool PointInsideRectangle(Point p, int numEl)
        {
            throw new NotImplementedException();
        }
                
        private double GetValueInPointInTriangle(Point p, int NumEl)
        {
            double detD;

            int[] _elem;

            _elem = _nvtr[NumEl];
            double a01, a02, a03, a11, a12, a13, a21, a22, a23;

            detD = (_points[_elem[1]].X - _points[_elem[0]].X) * (_points[_elem[2]].Y - _points[_elem[0]].Y) - (_points[_elem[2]].X - _points[_elem[0]].X) * (_points[_elem[1]].Y - _points[_elem[0]].Y);

            a01 = (_points[_elem[1]].X * _points[_elem[2]].Y - _points[_elem[2]].X * _points[_elem[1]].Y) / detD;
            a02 = (_points[_elem[2]].X * _points[_elem[0]].Y - _points[_elem[0]].X * _points[_elem[2]].Y) / detD;
            a03 = (_points[_elem[0]].X * _points[_elem[1]].Y - _points[_elem[1]].X * _points[_elem[0]].Y) / detD;

            a11 = (_points[_elem[1]].Y - _points[_elem[2]].Y) / detD;
            a12 = (_points[_elem[2]].Y - _points[_elem[0]].Y) / detD;
            a13 = (_points[_elem[0]].Y - _points[_elem[1]].Y) / detD;

            a21 = (_points[_elem[2]].X - _points[_elem[1]].X) / detD;
            a22 = (_points[_elem[0]].X - _points[_elem[2]].X) / detD;
            a23 = (_points[_elem[1]].X - _points[_elem[0]].X) / detD;

            return _q[_elem[0]] * (a01 + a11 * p.X + a21 * p.Y) +
                   _q[_elem[1]] * (a02 + a12 * p.X + a22 * p.Y) +
                   _q[_elem[2]] * (a03 + a13 * p.X + a23 * p.Y);
        }

        private double GetValueInPointInRectangle(Point p, int NumEl)
        {
            int[] _elem;
            _elem = _nvtr[NumEl];

            Point
                p0 = _points[_elem[0]],
                p1 = _points[_elem[1]],
                p2 = _points[_elem[2]];

            double hx = p1.X - p0.X;
            double hy = p2.Y - p0.Y;

            return  _q[_elem[0]] * ((p1.X - p.X)/hx) * ((p2.Y - p.Y)/hy) +
                    _q[_elem[1]] * ((p.X - p0.X)/hx) * ((p2.Y - p.Y)/hy) +
                    _q[_elem[2]] * ((p1.X - p.X)/hx) * ((p.Y - p0.Y)/hy) +
                    _q[_elem[3]] * ((p.X - p0.X)/hx) * ((p.Y - p0.Y)/hy);
        }

        private static bool PointInsideTriangle(Point p, Point p1, Point p2, Point p3)
        {
            double ab, ac, bc, ap, bp, cp;
            double pp, sABC;
            double s1, s2, s3;

            ab = Math.Sqrt(Math.Pow(p1.X - p2.X, 2) + Math.Pow(p1.Y - p2.Y, 2));
            ac = Math.Sqrt(Math.Pow(p1.X - p3.X, 2) + Math.Pow(p1.Y - p3.Y, 2));
            bc = Math.Sqrt(Math.Pow(p3.X - p2.X, 2) + Math.Pow(p3.Y - p2.Y, 2));
            ap = Math.Sqrt(Math.Pow(p1.X - p.X, 2) + Math.Pow(p1.Y - p.Y, 2));
            bp = Math.Sqrt(Math.Pow(p2.X - p.X, 2) + Math.Pow(p2.Y - p.Y, 2));
            cp = Math.Sqrt(Math.Pow(p3.X - p.X, 2) + Math.Pow(p3.Y - p.Y, 2));

            pp = (ab + ac + bc) / 2.0;
            sABC = Math.Sqrt(pp * (pp - ab) * (pp - ac) * (pp - bc));
            pp = (ab + ap + bp) / 2.0;
            s1 = Math.Sqrt(pp * (pp - ab) * (pp - ap) * (pp - bp));
            pp = (ac + ap + cp) / 2.0;
            s2 = Math.Sqrt(pp * (pp - ac) * (pp - ap) * (pp - cp));
            pp = (bc + bp + cp) / 2.0;
            s3 = Math.Sqrt(pp * (pp - bc) * (pp - bp) * (pp - cp));

            return Math.Abs(sABC - (s1 + s2 + s3)) < 1e-15;
        }

        private bool PointInsideTriangle(Point p, int NumEl)
        {
            int[] _elem = _nvtr[NumEl];
            return PointInsideTriangle(p, _points[_elem[0]], _points[_elem[1]], _points[_elem[2]]);            
        }

        private bool PointInsideRetangle(Point p, int NumEl)
        {
            int[] _elem = _nvtr[NumEl];
            double minX, maxX, minY, maxY;

            minX = minY = 10000D;
            maxX = maxY = -10000D;

            for (var i = 0; i < 4; ++i)
            {
                if (_points[_elem[i]].X < minX)
                    minX = _points[_elem[i]].X;

                if (_points[_elem[i]].X > maxX)
                    maxX = _points[_elem[i]].X;

                if (_points[_elem[i]].Y < minY)
                    minY = _points[_elem[i]].Y;

                if (_points[_elem[i]].Y > maxY)
                    maxY = _points[_elem[i]].Y;
            }

            if (p.X >= minX && p.X <= maxX && p.Y >= minY && p.Y <= maxY)
                return true;

            return false;
        }

        private bool PointInsideFiniteElement(Point p, int NumEl)
        {
            if (_elements == ElementType.Rectangle)
                return PointInsideTriangle(p, NumEl);
            else
                return PointInsideRetangle(p, NumEl);
        }

        /// <summary>
        /// Rebuild quad mesh
        /// </summary>
        private void RebuildMesh()
        {
            var nvtr = new List<int[]>();
            var q = _q.ToList();

            for(var i=0; i<_nvtr.Count; ++i)
            {
                var elem = _nvtr[i];
                var x = (_points[elem[1]].X + _points[elem[0]].X) / 2D;
                var y = (_points[elem[2]].Y + _points[elem[0]].Y) / 2D;

                var num = _points.Count;
                var p = new Point { X = x, Y = y };
                _points.Add(p);

                nvtr.Add(new[] { elem[0], elem[1], num, elem[4] });
                nvtr.Add(new[] { elem[0], num, elem[2], elem[4] });
                nvtr.Add(new[] { elem[0], elem[1], num, elem[4] });
                nvtr.Add(new[] { elem[2], elem[3], num, elem[4] });

                var s = GetValueInPointInRectangle(p, i);
                q.Add(s);
            }

            _q = q.ToArray();
            _nvtr = nvtr;
        }

        private void FragmentMesh()
        {
            var nvtr = new List<int[]>();
            var q = _q.ToList();

            for (var i = 0; i < _nvtr.Count; ++i)
            {
                var elem = _nvtr[i];
                var x = (_points[elem[0]].X + _points[elem[1]].X + _points[elem[2]].X) / 3D;
                var y = (_points[elem[0]].Y + _points[elem[1]].Y + _points[elem[2]].Y) / 3D;

                var num = _points.Count;
                var p = new Point { X = x, Y = y };
                _points.Add(p);

                nvtr.Add(new[] { elem[0], elem[1], num, elem[3] });
                nvtr.Add(new[] { elem[0], elem[2], num, elem[3] });
                nvtr.Add(new[] { elem[1], elem[2], num, elem[3] });

                var s = GetValueInPointInTriangle(p, i);
                q.Add(s);
            }

            _q = q.ToArray();
            _nvtr = nvtr;
        }

        VTreeNode HeadNode = new VTreeNode ();

        private void CreateVolumeTree (Volume mainVol, int deep)
        {
            HeadNode = new VTreeNode ();
            HeadNode.Vol = mainVol;
            
            
            if (deep > 0) {
                double x0 = mainVol.X0;
                double x1 = mainVol.X1;
                double y0 = mainVol.Y0;
                double y1 = mainVol.Y1;
                AddNodeToTree (HeadNode.Childs[0], new Volume ((x1 - x0) / 2, (y1 - y0) / 2, x1, y1), HeadNode, deep - 1);
                AddNodeToTree (HeadNode.Childs[1], new Volume ((x1 - x0) / 2, y0, x1, (y1 - y0) / 2), HeadNode, deep - 1);
                AddNodeToTree (HeadNode.Childs[2], new Volume (x0, (y1 - y0) / 2, (x1 - x0) / 2, y1), HeadNode, deep - 1);
                AddNodeToTree (HeadNode.Childs[3], new Volume (x0, y0, (x1 - x0) / 2, (y1 - y0) / 2), HeadNode, deep - 1);
            }
        }

        private void AddNodeToTree (VTreeNode node, Volume vol, VTreeNode parent, int deep)
        {
            node.Vol = vol;
            node.Parent = parent;

            if (deep > 0) {
                double x0 = vol.X0;
                double x1 = vol.X1;
                double y0 = vol.Y0;
                double y1 = vol.Y1;

                AddNodeToTree (HeadNode.Childs[0], new Volume ((x1 - x0) / 2, (y1 - y0) / 2, x1, y1), HeadNode, deep - 1);
                AddNodeToTree (HeadNode.Childs[1], new Volume ((x1 - x0) / 2, y0, x1, (y1 - y0) / 2), HeadNode, deep - 1);
                AddNodeToTree (HeadNode.Childs[2], new Volume (x0, (y1 - y0) / 2, (x1 - x0) / 2, y1), HeadNode, deep - 1);
                AddNodeToTree (HeadNode.Childs[3], new Volume (x0, y0, (x1 - x0) / 2, (y1 - y0) / 2), HeadNode, deep - 1);
            }
        }

        private void IncludeElementsInTree ()
        {
            for (int i = 0; i < _nvtr.Count; i++) {
                int[] elem = _nvtr[i];
                ElementInTree (elem, HeadNode);
            }
        }

        private void ElementInTree (int [] elem, VTreeNode node)
        {
            int numNodes = elem.Count<int> ();
            //выбираем узлы элемента
            List<Point> nodes = new List<Point> ();
            foreach (int ind in elem) {
                nodes.Add (_points[ind]);
            }

            double x0 = node.Vol.X0;
            double x1 = node.Vol.X1;
            double y0 = node.Vol.Y0;
            double y1 = node.Vol.Y1;

            int ru = 0;
            int rd = 0;
            int lu = 0;
            int ld = 0;

            foreach (Point n in nodes) {
                if (n.X > (x1 - x0) / 2) { //правая половина
                    if (n.Y > (y1 - y0) / 2) { //правая верхняя
                        ru++;
                    }
                    else { //правая нижняя
                        rd++;
                    }
                }
                else { //левая половина
                    if (n.Y > (y1 - y0) / 2) { //левая верхняя
                        lu++;
                    }
                    else {
                        ld++;
                    }
                }
            }

            //переходим к следующей подобласти
            if (ru >= 2) {
                ElementInTree (elem, node.Childs[0]);
            }
            else if (rd >= 2) {
                ElementInTree (elem, node.Childs[1]);
            }
            else if (lu >= 2) {
                ElementInTree (elem, node.Childs[2]);
            }
            else if (ld >= 2) {
                ElementInTree (elem, node.Childs[3]);
            }
        }
        
    }

    public enum BasisType
    {
        Linear = 2,
        Quad = 4,
        Cube = 8
    }

    public enum ElementType
    {
        Triangle = 3,
        Rectangle = 4
    }
}
